/*
   gpuGard.cpp

   Created by Miroslav Karpis on 1/12/11.
   Copyright 2011 Miroslav Karpis. All rights reserved.

   This file is part of molecEvol.

   molecEvol is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   Foobar is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with Foobar.  If not, see <http://www.gnu.org/licenses/>.
 */

#include "gpuGard.h"

gpuGard::gpuGard()
{
    //initialize constant GARD variables
    //TODO - simulation parameters needs to be loaded from file

    //GARD variables
    mean = -4.0;
    sigma = 4.0;
    Kf = pow(10.0, -2.0);
    Kb = pow(10.0, -4.0);
    Ng = 5;
    Nmax = Ng;
    gardPi = 1.0/Ng;

    //general constants
    pi = 3.1415926535;
    euler = 2.71828;
}

//------
void gpuGard::initGard(cl_uint nMolecTypes, float *Bij)
{
    cout << "Initiazing population with: " << nMolecTypes << endl;

    //generate random (based on log. normal) amount of molecules for numOfMolecTypes
    srand(time(NULL));

    clock_t start, end;
    start = clock();

    for(cl_uint i=0; i<nMolecTypes*nMolecTypes; i++)
    {
        Bij[i] = rand_logNormal(0.0, 1.0);
    }

    end = clock();
    /*
    cout << "Time required for execution: "
         << (float)(end-start)/CLOCKS_PER_SEC
         << " seconds." << "\n\n";
*/
}

//------
bool gpuGard::runGard(float *Bij)
{
    //calculate joion/leave rate
    float rate[Ng];

    //first we need to find current size of assembly
    float N = 0.0;
    for(cl_uint i=0;i<Ng;i++)
    {
        N = N + gardAssemb.molecules[i];
    }

    cout << "total size of assembly is:" << N << endl;

    //-------START - JOIN/LEAVE RATE

    //loop through all 100 molec. types and find the join/leave rate
    //N - is assembly current size
    //ni - molecular counts within the assembly
    for(cl_uint i = 0; i<Ng; i++)
    {
        float BijSum = 0.0;
        //get sum Bij value
        for(cl_uint j=0; j<Ng; j++)
        {
           // cout << "Bij:" << Bij[i*Ng+j] << "---" << i*Ng+j
           //         << "---" << gardAssemb.molecules[j]/N;
            BijSum = BijSum + (Bij[i*Ng+j]*(gardAssemb.molecules[j]/N));
        }

        rate[i] = (Kf*gardPi*N - Kb*gardAssemb.molecules[i])*(1+BijSum);
        cout << " rate: "<< rate[i] << "sum: " << BijSum << endl;
    }

    //-------END - CALCULATE JOIN/LEAVE RATE

    return true;
}

//------
float gpuGard::rand_logNormal(float mean, float variance) {
    static float n2 = 0.0;
    static int n2_cached = 0;
    if (!n2_cached) {
        float x, y, r;
        do {
            //TODO random seed needs to be moved to config file
            x = (2.0*(rand() % 10000))/10000 - 1;
            y = (2.0*(rand() % 10000))/10000 - 1;

            r = x*x + y*y;
        } while (r == 0.0 || r > 1.0);

        {
            float d = sqrt(-2.0*log(r)/r);
            float n1 = x*d;
            n2 = y*d;

            //here we apply log-normal distrib.
            float result = exp(n1*variance + mean);

            n2_cached = 1;
            return result;
        }
    } else {
        n2_cached = 0;
        //here we apply log-normal distrib.
        return exp(n2*variance + mean);
    }
}

